Disorder-Induced Critical Phenomena in Hysteresis: 
Numerical Scaling in Three and Higher Dimensions 
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We present numerical simulations of avalanches and critical 
phenomena associated with hysteresis loops, modeled using 
the zero-temperature random-field fsing model. We study 
the transition between smooth hysteresis loops and loops with 
a sharp jump in the magnetization, as the disorder in our 
model is decreased. In a large region near the critical point, we 
find scaling and critical phenomena, which are well described 
by the results of an e expansion about six dimensions. We 
present the results of simulations in 3, 4, and 5 dimensions, 
with systems with up to a billion spins (1000 3 ). 



I. INTRODUCTION 

The increased interest in real materials in condensed 
matter physics has brought disordered systems into the 
spotlight. Dirt changes the free energy landscape of a 
system, and can introduce metastable states with large 
energy barriers pL This can lead to extremely slow re- 
laxation towards the equilibrium state. On long length 
scales and practical time scales, a system driven by an 
external field will move from one metastable local free- 
energy minimum to the next. The equilibrium, global 
free energy minimum and the thermal fluctuations that 
drive the system toward it, are in this case irrelevant. 
The state of the system will instead depend on its his- 
tory. 

The motion from one local minima to the next is a col- 
lective process involving many local (magnetic) domains 
in a local region - an avalanche. In magnetic materi- 
als, as the external magnetic field H is changed contin- 
uously, these avalanches lead to the magnetic noise: the 
Barkhausen effect This effect can be picked up as 

voltage pulses in a coil surrounding the magnet. The dis- 
tribution of pulse (avalanche) sizes is found §-§ to follow 
a power law with a cutoff after a few decades, and was 
interpreted by some |(| to be an example of self-organized 
criticality (SOC) j?). (In SOC, a system organizes itself 
into a critical state without the need to tune an exter- 
nal parameter.) Other systems can exhibit avalanches 
as well. Several examples where disorder may play a 
part are: superconducting vortex line avalanches ||, re- 
sistance avalanches in superconducting films M , and cap- 



illary condensation of helium in Nuclepore pfl ] . 

The history dependence of the state of the system leads 
to hysteresis. Experiments with magnetic tapes JiT[ have 
shown that the shape of the hysteresis curve changes 
with the annealing temperature. The hysteresis curve 
goes from smooth to discontinuous as the annealing tem- 
perature is increased. This transition can be explained 
in terms of a plain old critical point with two tunable 
parameters: the annealing temperature and the exter- 
nal field. At the critical temperature and field, the cor- 
relation length diverges, and the distribution of pulse 
(avalanche) sizes follows a power law. 

We have argued earlier |2^] that the Barkhausen 
noise experiments can be quantitatively explained by a 
model (lj] with two tunable parameters (external field 
and disorder), which exhibits universal, non-equilibrium 
collective behavior. The model is athermal and incor- 
porates collective behavior through nearest neighbor in- 
teractions. The role of dirt or disorder, as we call it, 
is played by random fields. This paper presents the re- 
sults and conclusions of a large scale simulation of that 
model: the non-equilibrium zero-temperature Random 
Field Ising Model (RFIM), with a deterministic dynam- 
ics. The results compare very well to our e expansion 
|l3],[l4j , and to experiments in Barkhausen noise J25| . 

We should mention that there are other models for 
avalanches in disordered magnets. There is a large body 
of work on depinning transitions and the motion of the 
single interface ||J,[lJ|,[l6|] . In these models, avalanches 
occur only at the growing interface. Our model though, 
deals with many interacting interfaces: avalanches can 
grow anywhere in the system. Models of hysteresis simi- 
lar to ours exist [p3| , including ones with random bonds 
|l7| , ^8| and random anisotropics. 

This paper is a condensed version of an unpublished 
manuscript, available electronically Jl9|j , We focus here 
on the numerical results and scaling methods in dimen- 
sions three through six. Some of the other topics touched 
upon in the original manuscript are being published sep- 
arately. Our interpretation of the behavior in dimen- 
sion two has been substantially altered by further anal- 
ysis [ p0[ . A full description of the numerical method is 
available, including sample code and executables on the 
Web . For a full discussion of the behavior in mean 
field theory, and interesting behavior below the critical 
point in seven and nine dimensions, we refer the reader 
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to the electronic version of the original manuscript [^9| 
and to recent work on the Bethe lattice |22| . 



II. THE MODEL 

The model we use is the zero-temperature random- 
field Ising model |23|j2^ , ^2] , p5[ | , which we briefly review 
here. Magnetic domains are represented by spins Sj on a 
hypercubic lattice, which can take two values: Sj = ±1 . 
The spins interact ferromagnetically with their nearest 
neighbors with a strength Jy , and are exposed to a uni- 
form magnetic field H (which is directed along the spins) . 
Disorder is simulated by a random field hi, associated 
with each site of the lattice, which is given by a gaussian 
distribution function p{hi): 



P{hi) 



1 



2nR 



e 



(1) 



of width proportional to R, which we call the disorder. 
The Hamiltonian is then 



H = — JjjSjSj - /](H + hj)sj 



(2) 



<hj> 



For the analytic calculation, as well as the simulation, 
we have set the interaction between the spins to be inde- 
pendent of the spins and equal to one for nearest neigh- 
bors, Jij = J = 1 , and zero otherwise. We use periodic 
boundary conditions in the results of this paper; we've 
checked that the results near R c are unchanged when a 
slab of pre-flipped spins is introduced (fixed boundary 
conditions along two sides). 

The dynamics is deterministic, and is defined such that 
a spin Si will flip only when its local effective field : 



J^Sj + H + hi 



(3) 



changes sign. All the spins start pointing down (sj = — f 
for all i). As the field is adiabatically increased, a spin 
will flip. Due to the nearest neighbor interaction, a 
flipped spin will push a neighbor to flip, which in turn 
might push another neighbor, and so on, thereby gener- 
ating an avalanche of spin flips. During each avalanche, 
the external field is kept constant. For large disorders, 
the distribution of random fields is wide, and spins will 
tend to flip independently of each other. Only small 
avalanches will exist, and the magnetization curve will 
be smooth. On the other hand, a small disorder im- 
plies a narrow random field distribution which allows 
larger avalanches to occur. As the disorder is lowered, 
at the disorder R — R c and field H = H c , an infinite 
avalanche in the thermodynamic system will occur for 
the first time, and the magnetization curve will show a 
discontinuity. Near R c and H c , we find critical scaling 



behavior and avalanches of all sizes. Therefore, the sys- 
tem has two tunable parameters: the external field H 
and the disorder R. We found from the mean field cal- 
culation [^3|jT^] and the simulation that a discontinuity 
in the magnetization exists for disorders R < R c , at the 
field H C (R) > H C (R C ), but that only at (R C ,H C ), do we 
have critical behavior. For finite size systems of length L, 
the transition occurs at the disorder R^" (L) near which 
avalanches first begin to span the system in one of the d 
dimensions (spanning avalanches). The effective critical 
disorder R e J' (L) is larger than R c , and R%" (L) — > R c as 
L — > oo. 

The algorithm we use to simulate this model is de- 
scribed in a separate manuscript pl[ . For a simulation 
with N spins, the computer time scales as NlogN and 
the memory required for the simulation scales to one bit 
per spin (i.e., we do not store the random fields). 



III. SCALING 

We use data obtained from the simulation to find and 
describe the critical transition. We do so using scaling 
collapses, which we review briefly here. For example, 
the magnetization as a function of external field H is 
expected to have the form 

M(H,R)-M C (H C ,R C ) ~ \rf M ± (h'/\rf s ) (4) 

where M c is the critical magnetization (the magneti- 
zation at H c , for R = R c ), r = (R — R c )/R and 
h = (H — H c ) are the reduced disorder and reduced field 
respectively, 



h' = h + Br 



(5) 



is a (non-universal) rotation between the experimental 
control variables (r,h) and the scaling variables (r,h'), 
and M± is a universal scaling function (± refers to 
the sign of r).[] Scaling is expected asymptotically for 
small r and h — i.e., for H near H c and R near R c . 
The critical exponent gives the scaling for the mag- 
netization at the critical field H c [h = 0). If we plot 
\r\-P(M(H, R)-M C (H C , R c ) versus h/r^, we should ob- 
tain the curve Ai(x), independently of what disorder R 
we choose (so long as it is close to R c ): different ex- 
perimental and numerical data sets should collapse onto 



1 In the plots shown in this paper, we use r = (R — R c )/R, 
which we have found produces better collapses than using 
r = (R — Rc)/ Rc- The latter is more traditional, but the two 
definitions agree as R — > R c , and differ by an amount which 
is irrelevant in a renormalization-group sense. One method 
we use to estimate error in our exponents is to compare ex- 
trapolations based on the two definitions. 
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one universal curve M. (x) . (Actually, one has two curves 
M± depending on whether R > R c or R < R c .) We use 
scaling forms similar to (0) to analyze all of our measure- 
ments. 

One can easily show using the scaling form (^) that the 
magnetization scales with a power law M — M c <~ h s at 
R c , and that the jump in the magnetization (the size of 



as one varies 



the infinite avalanche) scales as AM 
the disorder below R c . Thus the critical exponents j3 
and 5 give the power laws for the singularities in these 
measured quantities: indeed, that is how these exponents 
were originally defined and measured. In our system, 
we will find that directly measuring power laws is not 
effective in getting good exponents: the critical regime 
is so large that we need both to use the general scaling 
form and to extrapolate to the critical point. 

The explanatory power of the theory resides in the 
fact that the same universal critical exponents (3 and 
6 and the same universal function M.{x) should be ob- 
tained by simulations at different values of the disorder, 
simulations of different Hamiltonians, and simulations 
of real experiments, so long as the systems share cer- 
tain important features and symmetries (so long as they 
lie in the same universality class). The underlying ex- 
planation for why universality and scaling should occur 
near the critical point is given by the renormalization 
group [ f^ , pli||l4| , |l9| ] . Above six dimensions, fluctuations 
are asymptotically not important, and we can calculate 
Ai(x) and the values of and 6 from mean field the- 
ory (Pmf — 1/2, Smf = 3 f0| . Below six dimensions, 
the exponents and scaling curves are non-trivial, and to 
find them one must rely on either perturbative meth- 
ods plO], experiments, or numerical methods [ ^2|j25| 
as used here. 



IV. THE SIMULATION RESULTS 

The following measurements were obtained from the 
simulation as a function of disorder R: 

• the magnetization M(H, R) as a function of the 
external field H. 

• the avalanche size distribution integrated over the 
field H: D int (S,R). 

• the avalanche correlation function integrated over the 
field H: G int (x,R). 

• the number of spanning avalanches N(L, R) as a 
function of the system length L, integrated over the 
field H. 

• the discontinuity in the magnetization AM(L, R) as 
a function of the system length L. 

• the second (S 2 ) mt (L, R), third (S 3 ) mt (L, R), and 
fourth (S )i n t{L,R) moments of the avalanche size 
distribution as a function of the system length L, 
integrated over the field H. 

In addition, we have measured: 



• the avalanche size distribution D(S, H, R) as a 
function of the field H and disorder R. 

• the distribution of avalanche times D^^ (S, t) as a 
function of the avalanche size S, at R = R c , integrated 
over the field H . 



A. Magnetization Curves 

Unfortunately the most obvious measured quantity in 
our simulations, the magnetization curve M(£f), is the 
one which collapses least well in our simulations. We 
start with it nonetheless. 
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FIG. 1. (a) Magnetization curves in 3 dimensions for 

size L = 320, and three values of disorder. The curves are av- 
erages of up to 48 different random field configurations. Note 
the discontinuity in the magnetization for R = 2.20. In finite 
size systems, the discontinuity in the magnetization curve oc- 
curs even for R> R c (R c — 2.16 in 3 dimensions), (b) Scaling 
collapse (see text) of the magnetization curves in 3 dimen- 
sions for size L — 320. The disorders range from R — 2.35 to 
R — 3.20. The critical magnetization is chosen as M c — 0.9 
from an analysis of the magnetization curves and is kept fixed 
during the collapse. The universal exponents are f3 = 0.036, 
135 — 1.81. The non-universal critical field H c = 1.435, criti- 
cal disorder R c = 2.16, and rotation parameter B — 0.39. 



3 



Figure [l]a shows the magnetization curves obtained 
from our simulation in 3 dimensions for several values 
of the disorder R. As the disorder R is decreased, a dis- 
continuity or jump in the magnetization curve appears 
where a single avalanche occupies a large fraction of the 
total system. In the thermodynamic limit this would be 
the infinite avalanche: the largest disorder at which it 
occurs is the critical disorder R c . For finite size systems, 
like the ones we use in our simulation, we observe an 
avalanche which spans the system at a higher disorder, 
which gradually approaches R c as the system size grows. 
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FIG. 2. dM/dH curves in 3 dimensions (a) Derivative of 
the magnetization M with respect to the field H for disorders 
R = 2.35, 2.4, 2.45, 2.5, 2.6, 2.7, 2.85, 3.0, and 3.2 (highest 
to lowest peak), (b) Scaling collapse of the data in (a) with 
P = 0.036, p6 = 1.81, B = 0.39, H c = 1.435, and R c = 2.16. 
While the curves are not collapsing onto a single curve, the 
quality of the collapse is quite similar to that found at similar 
distances from R c in mean field theory Jis| , for which we know 
analytically that scaling works as R — » R c . 

Figure ^ shows the slope dM/dH and its scaling col- 
lapse. By using this derivative, the critical region is em- 
phasized as the peak in the curve, and the dependence 
on the parameter M c drops out. The lower graphs in 
figures |l|b and ||b show the scaling collapses of the mag- 



netization and its slope. Clearly in neither case is all the 
data collapsing onto a single curve. This would be dis- 
tressing, were it not for the fact that this also occurs in 
mean field theory |l9| at a similar distance to the critical 
point. 

Because the scaling of the magnetization is so bad, we 
use other quantities to estimate the critical exponents 
and the location of the critical point (tables | and p| ). 
Fixing these quantities, we use the collapse of the dM / dH 
curves to extract the rotation B mixing the experimental 
variables r and h into the scaling variable h' — h + Br 
(equation ||) . 



B. Avalanche Size Distribution 



1. Integrated Avalanche Size Distribution 
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FIG. 3. Integrated avalanche size distribution 
curves in 3 dimensions for 320 3 spins and disorders 
R = 4.0,3.2, and 2.6. The last curve is at R = 2.25, for 
a 1000 3 spin system. The 320 3 curves are averages over up 
to 16 initial random field configurations. The inset shows the 
scaling collapse of the integrated avalanche size distribution 
curves in 3 dimensions, using r = (R — R c ) / R, r + apS = 2.03, 
and a = 0.24, for sizes 160 3 , 320 3 , 800 3 , and 1000 3 , and disor- 
ders ranging from R = 2.25 to R = 3.2 (Rc = 2.16). The two 
top curves in the collapse, at R — 3.2, show noticeable correc- 
tions to scaling. The thick dark curve through the collapse is 
the fit to the data (see text). In the main figure, the distri- 
bution curves obtained from the fit to the collapsed data are 
plotted (thin lines) alongside the raw data (thick lines). The 
straight dashed line is the expected asymptotic power law be- 
havior: S -2 ' 03 , which does not agree with the measured slope 
of the raw data due to the shape of the scaling function (see 
text). 

In our model the spins flip in avalanches: each spin 
can kick over one or more neighbors in a cascade. 
These avalanches come in different sizes. The integrated 
avalanche size distribution is the size distribution of all 
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the avalanches that occur in one branch of the hystere- 
sis loop (for H from — oo to oo). Figure || (25) shows 
some of the raw data (thick lines) in 3 dimensions. Note 
that the curves follow an approximate power law behav- 
ior over several decades. Even 50% away from criticality 
(at R = 3.2), there are still two decades of scaling, which 
implies that the critical region is large. In experiments, 
a few decades of scaling could be interpreted in terms of 
self-organized criticality (SOC). However, our model and 
simulation suggest that several decades of power law scal- 
ing can still be present rather far from the critical point 
(note that the size of the critical region is non-universal). 
The slope of the log-log avalanche size distribution at R c 
gives the critical exponent r+a(3S. Notice, however, that 
the apparent slopes in figure || continue to change even 
after several decades of apparent scaling is obtained. The 
cutoff in the power law diverges as the critical disorder 
R c is approached. This cutoff size scales as S ~ |r| _1 / CT . 

These critical exponents can be obtained by using a 
scaling collapse for the curves of figure ^, shown in the 
inset. The scaling form is 

D mt (S,R) ~ S- {T+af) ^ V { l nt) {S°\r\) (6) 

where r D^ nt ' is the scaling function (the + sign indicates 
that the collapsed curves are for R > R c ). We are suffi- 
ciently far from the critical point that corrections to scal- 
ing are important: as described in reference ]l9|] , we do 
collapses for small ranges of R and then linearly extrap- 
olate the best-fit critical exponents to R c . We estimate 
from this curve that the critical exponents t + <j(38 = 2.03 
and a = 0.24 

The scaling function V% nt) (X) with X = S a \r\ is a 
universal prediction of our model. To facilitate compar- 
isons with experiments, we fit a curve to the data collapse 
in the inset of figure |3|. We have fit the scaling collapses 
in dimensions 3, 4, and 5 to a phenomenological form of 
an exponential times a polynomial. In three dimensions, 
our fit is 

V { l nt) {X) = e -°- n9xl " x 
(0.021 + 0.002X + 0.53LX 2 - 0.266X 3 + 0.261X 4 ) (7) 

where 1/er = 4.20. The distribution curves obtained us- 
ing the above fit are plotted (thin lines in figure ||) along- 
side the raw data (thick lines). They agree remarkably 
well even far above R c . We should recall though, that 
the fitted curve to the collapsed data can differ from the 
"real" scaling function even for large sizes and close to 
the critical disorder (in mean field |l9f| the error in the 
corresponding curve was about 10%). 

The scaling function in the inset of figure || has a pe- 
culiar shape: it grows by a factor of ten before cutting 
off. The consequence of this bump in the shape is that in 
the simulations it takes many decades in the size distri- 
bution for the slope to converge to the asymptotic power 



law. This can be seen from the comparison between a 
straight line fit through the R = 2.25 (billion spin) sim- 
ulation in figure ||| and the asymptotic power law S~ 2m 
obtained from extrapolating the scaling collapses (thick 
dashed straight line in the same figure). A similar bump 
exists in other dimensions and mean field as well. Fig- 
ure U shows the scaling functions in different dimensions 
and in mean field. In this graph, the scaling functions are 
normalized to one and the peaks are aligned (the scaling 
forms allow this). The curves plotted in figure |] are not 
raw data but fits to the scaling collapse in each dimen- 
sions, as was done in the inset of figure [| For 5,4, and 
2 dimensions, we have respectively: 

V { 5 mt) (X) = e -°- 51&xl " x 
(0.112 + 0.459X - 0.260X 2 + 0.201X 3 - 0.050X 4 ) (8) 

vi mt) (X) = e -°- 954xl/ " x 
(0.058 + 0.396X + 0.248X 2 - 0.140X 3 + 0.026X 4 ) (9) 

V { 2 mt) (X) = e-l™ 17 ' X 

(0.492 - 4A72X + 14.702X 2 - 

20.936X 3 + 11.303X 4 ) (10) 

with l/a = 2.35,3.20, and 10.0. The errors in the fits 
are again in the 10% range, judging from mean-field the- 
ory QJ. 




S°lrl 

FIG. 4. Integrated avalanche size distribution scal- 
ing functions in 2, 3, 4, and 5 dimensions, and mean 
field. The curves are fits (see text) to the scaling collapses 
done with exponents from Table | and corresponding calcula- 
tions in two dimensions |20| . The peaks are aligned to fall on 
(1,1). Due to the "bump" in the scaling function the power 
law exponent can not be extracted from a linear fit to the raw 
data for reasonable simulation sizes [E5|. 
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In mean field theory (dimensions six and greater) a 
similar fit |ll| to the analytical form of the scaling func- 
tion above R c gives 



V 



(int) 
MF 



(0.204 + 0.482X - 0.391X 2 
0.204X 3 - 0.048X 4 ) 



(11) 



It is clear from the figure that the growing bump in the 
scaling curves as the dimension decreases is a foreshad- 
owing of a zero in the scaling curve in two dimensions: 
this will be discussed further in reference |E0l. 



2. Binned in H Avalanche Size Distribution 

The avalanche size distribution can also be measured at 
a field H or in a small range of fields centered around H . 
We have measured this binned in H avalanche size distri- 
bution for systems at the critical disorder R c (r = 0). To 
obtain the scaling form, we start from the distribution of 
avalanches at field H and disorder R 



D(S,R,H) ~ S~ T V±(S° 



\h\/\r 



(12) 



where as before T>± is the scaling function and ± indi- 
cates the sign of r. 

The parameter B of equation |5|, which rotates the mea- 
sured axes (r, h) into the scaling axes (r, h' = h + Br), 
will be important (ll| only for large avalanches of size 
S > hr 1 ! near the critical point. In three and four di- 
mensions, this does not affect our scaling collapses; in 
five dimensions we account for it Jl9{ . 

The scaling function can be rewritten as 

T> ± (s cr \r\,(S' T \r\f 5 \h\/\rf 6 ^, where D± is a new scal- 
ing function and ± represents whether H is greater than 
or less than H c (i.e., H rather than R). Letting R — > R c , 
the scaling for the avalanche size distribution at the field 
H , measured at the critical disorder R r is: 



D(S,H)~S- T V ± (\h\S a0S ) 



(13) 



Figure |^a shows the binned in H avalanche size dis- 
tribution curves in 4 dimensions, for values of H below 
the critical field H c . (The curves and analysis are similar 
in 3 and 5 dimensions; results in 4 dimensions are used 
here for variety.) The simulation was done at the best 
estimate of the critical disorder R c (4.1 in 4 dimensions). 
The binning in H is logarithmic and started from an 
approximate critical field H c obtained from the magne- 
tization curves; better estimates of H c are then obtained 
from the binned distribution data curves and their col- 
lapses. Our best estimate for the critical field H c in 4 
dimensions is 1.265 ±0.007. The scaling form for the log- 
arithmically binned data is the same as in equation (13), 
if the log-binned data is normalized by the size of the 
bin. Figure ||b shows the scaling collapse for our data, 



both below and above the critical field H c . The "top" col- 
lapse gives the shape of the £>_ (H < H c ) function, while 
the "bottom" collapse gives the T>+ (H > H c ) function. 
Above the critical field H c , there are spanning avalanches 
in the system ]2§| ] . These are not included in the binned 
avalanche size distribution collapse shown in figure |B|b. 
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FIG. 5. (a) Binned in H avalanche size distribution 
in 4 dimensions for a system of 80 4 spins at R = 4.09 
(R c = 4.10). The critical field is H c = 1.265. The curves are 
averages over close to 60 random field configuration. Only 
a few curves are shown, (b) Scaling collapse of the binned 
avalanche size distribution for H < H c (upper collapse) and 
H > H c (lower collapse). The critical exponents are r = 1.53 
and a(3S = 0.54, and the critical field is H c = 1.265. The bins 
are at fields 1.162, 1.185, 1.204, 1.220, 1.234, 1.245, 1.254, 
1.276, 1.285, 1.296, 1.310, 1.326, 1.345, and 1.368. 

The exponent r which gives the power law behavior of 
the binned avalanche size distribution is obtained from 
collapses of neighboring curves as described above Jl9| , 
extrapolating to H — H c . The exponent a (38 is found to 
be very sensitive to H c , while r is not. We have there- 
fore used the values of r + <j(38 and a from the inte- 
grated avalanche size distribution collapses, and r from 
the binned avalanche size distribution collapses to further 
constrain H c (by constraining a (38), and to calculate (38. 



G 



The latter is then used to obtain collapses of the mag- 
netization curves. We should mention here that H c in 
all the dimensions is difficult to find and that it is influ- 
enced by finite sizes. The values listed in Table III are 



C. Avalanche Correlations 



the best estimates obtained from the largest system sizes 
we have. Nevertheless, systematic errors for H c could be 
larger than the errors given in Table [II. These errors 



could produce systematic errors for afld which depends 
on H c , and for (38 which is calculated from afiS: hence 
errors in these exponents could also be larger than the 
errors listed in Table 0. 
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FIG. 6. Linear fit to binned avalanche size distribu- 
tion curve in 4 dimensions, for a system of 80 4 spins at 
R c = 4.09. The magnetic field is H = 1.265. The straight 
solid line is a linear fit to the data for S < 13, 000 spins. The 
slope from the fit is 1.55 (this varies by not more than 3% 
as the range over which the fit is done is changed), while the 
exponent r obtained from the collapses and the extrapolation 
in figure | is 1.53 ± 0.08. 

From figure |^b, we see that the two binned avalanche 
size distribution scaling function do not have a "bump" as 
does the scaling function for the integrated avalanche size 
distribution (inset in figure [}]). Therefore, we expect that 
the exponent r which gives the slope of the distribution 
in figure ^a can also be obtained by a linear fit through 
the data curve closest to the critical field. Figure |^ shows 
the curve for the H — 1.265 bin (dashed curve) as well 
as the linear fit. The slope from the linear fit is 1.55 
while the value of r obtained from the collapses and the 
extrapolation in figure [| is 1.53 ± 0.08. 
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FIG. 7. Avalanche correlation function integrated 
over the field H in 3 dimensions, for L — 320. The 
curves are averages of up to 19 random field configurations. 
The critical disorder R c is 2.16. 

The avalanche correlation function G(x, R, H) mea- 
sures the probability that the initial spin of an avalanche 
will trigger, in that avalanche, another spin a distance 
x away. From the renormalization group description 
p3| , ^4[ , close to the critical point and for large distances 
x, the correlation function is given by: 



G(x,R, H) 



1 



d-2+7) 



G±(x/ar, h)) 



(14) 



where r and h are respectively the reduced disorder and 
field, G± (± indicates the sign of r) is the scaling function, 
d is the dimension, £ is the correlation length, and rj is 
called the "anomalous dimension" . Corrections can be 
shown to be subdominant p9| . The correlation length 
£ (r, h) is a macroscopic length scale in the system which 
is on the order of the mean linear extent of the largest 
avalanches. At the critical field H c (h=0) and near R c , 
the correlation length scales like £ ~ \ f\~ v > while for small 
field h it is given by 



t~\r\- v y ± (h/\rf s ) 



(15) 



where y± is a universal scaling function. The avalanche 
correlation function should not be confused with the clus- 
ter or "spin-spin" correlation which measures the prob- 
ability that two spins a distance x away have the same 
value. (The algebraic decay for this other, spin-spin cor- 
relation function at the critical point (r = and h = 0), 
is l/x d - i+f > Q.) 

We've mostly used, for historical reasons, a slightly 
different avalanche correlation function, which scales the 
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same way as the "triggered" correlation function G de- 
scribed above. Our function basically ignores the differ- 
ence between the triggering spin and the other spins in 
the avalanche: alternatively it calculates for avalanches 
of size S the correlation function for pairs of spins, 
and then averages over all avalanches (weighting each 
avalanche equally). We've checked that our correlation 
function agrees to within 3% with the "triggered" cor- 
relation function described above, for R > R c in three 
dimensions and above. (In two dimensions, the two def- 
initions differ more substantially but appear to scale in 
the same way pC|.) 

We have measured the avalanche correlation function 
integrated over the field H, for R > R c . For every 
avalanche that occurs between H — — oo and H = ±00, 
we keep a count on the number of times a distance x 
occurs in the avalanche. To decrease the computational 
time not every pair of spins is selected; instead we do 
a statistical sampling ]2lf| . The spanning avalanches are 
not included in our correlation measurement. Figure Q 
shows several avalanche correlation curves in 3 dimen- 
sions for L — 320. The scaling form for the avalanche 
correlation function integrated over the field H, close to 
the critical point and for large distances x, is obtained 
by integrating equation (|l4|): 

G mt (x,R) ~ J Q±(x/£(r,hj) dh (16) 

Using equation ( |l5j ) and defining u = h/\r\^ s , equation 
( |l6| ) becomes: 

G M (x,R)~ \rf s x-( d - 2+ ^ J Q±(x/\r\- v y±{u)) du 

(17) 

The integral (I) in equation ( |l7j ) is a function of x\r\ v 
and can be written as: 

1= {x\r\")-W" G ± (x\r\ v ) (18) 

to obtain the scaling form: 

G int{x,R)~^j;g ± {x\r\ v ) (19) 

where we have used the scaling relation (2 — if)v = (38 — (3 
(see Q for the derivation). 
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FIG. 8. Scaling collapse of the avalanche correlation 
function integrated over the field H, in 3 dimensions 

for L — 320. The values of the disorders range from R = 2.35 
to R — 3.0, with R c — 2.16. The exponents used in the 
collapse are v = 1.39 ±0.20 and d + P/v = 3.07 ±0.30. When 
collapses of neighboring curves are extrapolated to R c , we get 
a slightly smaller value of v = 1.37 ± 0.18. 
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FIG. 9. Anisotropics in the avalanche correlation 
function. The curves are for a system of 320 3 spins at 
R — 2.35. Four curves are shown on the graph: one is the 
avalanche correlation function integrated over the field H (as 
in figure [j]) , while the other three are measurements of the cor- 
relation along the three axis, the six face diagonals, and the 
four body diagonals. Avalanches involving more than four 
spins show no noticeable anisotropy: the critical point ap- 
pears to have spherical symmetry. The same result is found 
in 2 dimensions. 



Figure g shows the integrated avalanche correlation 
curves collapse in 3 dimensions for L = 320 and R > 
R c . The exponent v is obtained from such collapses 
by extrapolating to R = R c as was done for other col- 
lapses [^9) . The exponent (3/v can be obtained from these 
collapses too, but it is much better estimated from the 
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magnetization discontinuity covered below. The value 
of P/v listed in Table | is derived exclusively from the 
magnetization discontinuity collapses. 

We have also looked for possible anisotropies in the 
integrated avalanche correlation function in 2 and 3 di- 
mensions. The anisotropic integrated avalanche corre- 
lation functions are measured along "generalized diago- 
nals": one along the three axis, the second along the six 
face diagonals, and the third along the four body diag- 
onals. We compare the integrated avalanche correlation 
function and the anisotropic integrated avalanche corre- 
lation functions to each other, and find no anisotropics 
in the correlation, as can be seen from figure |p. 



D. Spanning Avalanches 
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FIG. 10. Spanning avalanches in 3, 4, and 5 dimen- 
sions, (a) Number of spanning avalanches TV in 3 di- 
mensions, occurring in the system between H = — oo to 
H — oo, as a function of the disorder R, for linear sizes L: 
20 (dot-dashed), 40 (long dashed), 80 (dashed), 160 (dotted), 
and 320 (solid). The critical disorder R c is at 2.16. The er- 
ror bars for each curve tend to be smaller than the error bar 
shown at the peak for disorders above the peak and larger for 
disorders below the peak. They are not given here for clarity. 
Note that the number of avalanches increases only slightly as 
the size is increased, (b) Number of spanning avalanches 
in 4 dimensions. The critical disorder is 4.1. (c) Number 
of spanning avalanches in 5 dimensions. The critical 
disorder is 5.96. Both in 4 and 5 dimensions, the peaks grow 
and shift towards R c as the size of the system is increased, (d) 
Collapse of the spanning avalanche curves in 4 dimen- 
sions for linear sizes L = 20, 40, and 80. The exponents are 
6 = 0.32 and v — 0.89, and the critical disorder is R c = 4.10. 
The collapse is done using r = (R c — R)/R. 

The critical disorder R c was defined earlier as the dis- 
order R at which an infinite avalanche first appears in 
the system, in the thermodynamic limit, as the disor- 
der is lowered. At that point, the magnetization curve 
will show a discontinuity at the magnetization M C (R C ) 
and field H C (R C ). For each disorder R below the critical 
disorder, there is one infinite avalanche that occurs at 
a critical field H C (R) > H C {R C ) while above R c 

there are only finite avalanches. This is the behavior for 
an infinite size system. In a finite size system far below 
and above R c the above picture is still true, but close to 
the critical disorder, as we approach the transition, the 
avalanches get larger and larger, and there will be a first 
point where one of them will span the system from one 
side to another in at least one direction. This avalanche 
is not the infinite avalanche; if the system was larger, 
this avalanche would typically be non-system spanning. 
Such an avalanche (which spans the system) we call a 
spanning avalanche. 

fn our numerical simulation, we find that for finite sizes 
L, there are not one but many such avalanches in 4 and 
5 dimensions (and maybe 3), and that their number in- 
creases as the system size increases p9[ |. Figures flO|(a-c) 
show the number of spanning avalanches as a function of 
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disorder _R, for different sizes and dimensions. In 4 and 5 
dimensions, the spanning avalanche curves become more 
narrow as the system size is increased. Also, the peaks 
shift toward the critical value of the disorder (4.1 and 5.96 
respectively), and the number of spanning avalanches at 
R c increases. This suggests that in 4 and 5 dimensions, 
for L — ► oo, there will be one infinite avalanche below R Cl 
none above, and an infinite number of infinite, spanning 
avalanches at the critical disorder R c . In 3 dimensions, 
the results are not conclusive, as noted both from fig- 
ure 10 1 and from the value of the spanning avalanche 
exponent 9 = 0.15 ± 0.15 defined below: a value of 9 = 
is consistent with one infinite or spanning avalanche at 
R c as L — ► oo. It is clear that 9 = in two dimensions, 
since spanning avalanches can't interpenetrate: it's thus 
plausible that 9 is near zero in three dimensions because 
it must vanish one dimension lower. 

In percolation, a similar multiplicity of infinite clus- 
ters f|(],|3l| as the system size is increased is found for 
dimensions above six which is the upper critical dimen- 
sion (UCD). The UCD is the dimension at and above 
which the mean field exponents are valid. Below six di- 
mensions, there is only one such infinite cluster in per- 
colation. The existence of a diverging number of infinite 
clusters in percolation is associated with the breakdown 
of the hyperscaling relation above six dimensions. Since 
a hyperscaling relation is a relation between critical ex- 
ponents that includes the dimension d of the system, it 
is always only satisfied up to and including the upper 
critical dimension. In our system, the upper critical di- 
mension is also six, but we find spanning avalanches in 
dimensions even below that. In a comment by Maritan et 
al. p9| ], it was suggested that our system should satisfy 
the hyperscaling relation dv — f3 = 1 /a found in percola- 
tion pip . But since our system has spanning avalanches 
below the upper critical dimension, this hyperscaling re- 
lation breaks down below six dimensions. Due to the 
existence of many spanning avalanches near R c , the new 
"violation of hyperscaling" relation for dimensions three 
and above becomes [fL3f: 



(d-6)v-/3= 1/a 



(20) 



where 9 is the "breakdown of hyperscaling" or spanning 
avalanches exponent defined below. One can check that 
our exponents in 3, 4, and 5 dimensions and mean field 
satisfy this equation (see Tables Q and |J). 

For the simulation, we define a spanning avalanche to 
be an avalanche that spans the system in a particular 
direction. We average over all the directions to obtain 
better statistics. Depending on the size and dimension 
of the system and the distance from the critical disorder, 
the number of spanning avalanches for a particular value 
of disorder R is obtained by averaging over as few as 5 
to as many as 2000 different random field configurations. 
We define the exponent 9 such that the number N of 



spanning avalanches, at the critical disorder R c , increases 
with the linear system size as: N ~ L e (9 > 0). The 
finite size scaling form p^j ] for the number of spanning 
avalanches close to the critical disorder is: 



N(L,R) ~ L e M ± {L 1 / V \r\) 



(21) 



where v is the correlation length exponent and Af± is the 
corresponding scaling function (± indicates the sign of 
r). The collapse is shown in figure [Tc|d. (We show the 
collapses in 4 dimensions here since the existence of span- 
ning avalanches in 3 dimensions is not conclusive.) These 
values are used along with the results from other collapses 
to obtain Table |. In the analysis of the avalanche size 
distribution, magnetization, and correlation functions for 
R > R c , how close we chose to come to the critical disor- 
der R c was determined by the spanning avalanches: we 
include no values R below the first value which exhibited 
a spanning avalanche. 



E. Magnetization Discontinuity 

We have mentioned earlier that in the thermodynamic 
limit, at and below the critical disorder R c , there is 
a critical field H C (R) > H C (R C ) at which the infinite 
avalanche occurs. Close to the critical transition, for 
small r < 0, the change in the magnetization due to 
the infinite avalanche scales as (equation (El)): 



AM (R) 



,/3 



(22) 



where r = (R c — R)/R, while above the transition, there 
is no infinite avalanche. 

In finite size systems, the transition is not as sharp: 
we have spanning avalanches above the critical disorder. 
If we measure the change in the magnetization due to all 
the spanning avalanches as a function of disorder R at 
various system sizes L, we expect it will obey finite-size 
scaling (as did the number of spanning avalanches): 



AM(L,R) ~ \rf AM±{L 1 ^\ 



(23) 



where AA4± is a universal scaling function. (The param- 
eter B here (equation ^|) is unimportant [fl9f because we 
AM is measured at h! — 0.) Defining a new universal 
scaling function AA4±: 

AM±(L 1 ^\r\) = {L x / V \r\)- AM±(L 1 ^ u \r\) (24) 

we obtain the scaling form: 

AM(L,R) ~ L~®l v AM±{L 1 ^\r\) (25) 
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F. Moments of the Avalanche Size Distribution 
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FIG. 11. Jump in the magnetization, in 4 dimen- 
sions (a) Change in the magnetization due to the spanning 
avalanches in 4 dimensions, for several linear sizes L, as a 
function of the disorder R. (b) Scaling collapse of the curves 
in (a) using r — (R c — R)/R. The exponents are \/v — 1.12 
and P/v = 0.19, and the critical disorder is R c = 4.1. 



Figures 11a and ^TJd show the change in the magneti- 
zation due to the spanning avalanches in 4 dimensions, 
and a scaling collapse of that data (similar results exist 
in 3 and 5 dimensions). Notice that as the system size 
increases, the curves approach the \r\" behavior. The ex- 
ponents 1/v and (3/u are extracted from scaling collapses 
(figure p"l"|b) and extrapolated to R c [|l9| . The value of (3 
is calculated from (3/v and the knowledge of v, and is 
the value used for collapses of the magnetization curves 
(discussed earlier). 
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FIG. 12. Second moments, (a) Second moments of the 
avalanche size distribution integrated over the field H, in 5 di- 
mensions. Error bars are largest for smaller disorders (shown 
on the curves). The curves have between 24 and 50 points, 
and the value of the second moment for each disorder is av- 
eraged over 3 to 100 different random field configurations, 
(b) Scaling collapse of the L = 10, 20, and 30 curves from (a) 
using r — (R c — R)/R. The exponents are 1/u = 1.47 and 
p = — (t + af3S — 3) /(tv = 2.95, and the critical disorder is 
R c = 5.96. 

The second moment of the integrated avalanche size 
distribution has a finite-size scaling form 



(26) 



where L is the linear size of the system, r is the re- 

— (2) 

duced disorder, 5^ is the scaling function, and v is 
the correlation length exponent. We can similarly de- 
fine the third and fourth moment, with the exponent 

— (t + cr(35 — 3)/au replaced by — (r + a (35 — 4)/av and 

— (r + a [35 — h)/av respectively. Figures [l2|a and p~2| h> 
show the second moments data in 5 dimensions for sizes 
L = 5, 10, 20, and 30, and a collapse (again, results in 3 
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and 4 dimensions are similar and we have chosen to show 
the curves in 5 dimensions for variety). The curves are 
normalized by the average avalanche size integrated over 
all fields H: S D(S, R, H, L) dS dH. The span- 

ning avalanches are not included in the calculation of the 
moments. We omit the L = 5 curve from the collapse; 
it doesn't collapse with the others well, presumably be- 
cause of subdominant finite size effects. The exponents 
for the third and fourth moment can be calculated from 
those of the second moment, and we find that they agree 
with the values obtained from their respective collapses. 



G. Avalanche Time Measurement 

The exponents we have measured so far are static scal- 
ing exponents: they do not depend on the dynamics of 
the model. If we measure the time an avalanche takes 
to occur, we are making a dynamical measurement. The 
time measurement in the numerical simulation is done by 
increasing the time clock by one for each shell of spins 
in the avalanche. That is, we implement time as a syn- 
chronous dynamics, where in each time step all unstable 
spins from the previous step are flipped. The scaling re- 
lation between the time t it takes an avalanche to occur 
and the linear size £ of the avalanche defines the critical 
exponent z p2l.p3|: 



(27) 



The exponent z is known as the dynamical critical ex- 
ponent. Equation ( B?j ) gives the scaling for the time it 
takes for a spin to "feel" the effect of another a distance £ 
away. Since the correlation length £ scales like r~ v close 
to the critical disorder, and the characteristic size S as 
r -1 / ", the time t then scales with avalanche size as: 



t ~ S a 



(28) 



In our simulation, we measure the distribution of times 
for each avalanche size S. The distribution of times 
D t (S, R, H, t) for an avalanche of size S close to the crit- 
ical field H c and critical disorder R r is 



D t (S, R, H,t) ~ S 



- q v M( S o 



\r\,h/\rf S ,t/S™ z ) (29) 



where q = r + cruz, and is denned such that 

/+oo pOO 
/ D t (S,R,H,t) dH dt = 
oo Jl 



(30) 



- 1 

where T>± was defined in the integrated avalanche size 
distribution section. The avalanche time distribution in- 
tegrated over the field H, at the critical disorder (r = 0) 
is: 

Df nt) {S,t) - rlT+^+H/™^ 1 ")^") (3i) 
which obtained from equation (E9j) (reference fl8|]). 




FIG. 13. (a) Avalanche time distribution curves in 
3 dimensions, for avalanche size bins from about 2000 to 
40000 spins (from upper left to lower right corner). The sys- 
tem size is 800 3 at R = 2.26. The curves are from only one 
random field configuration, (b) Scaling collapse of curves 
in (a). The values of the exponents are avz — 0.57 and 
(r + a/35 + avz) javz = 4.0. 



Figures |13Ja and |13| b show the avalanche time distri- 
bution integrated over the field H for different avalanche 
sizes, and a collapse of these curves using the above scal- 
ing form, for a 800 3 system at R = 2.260 (just above 
the range where spanning avalanches occur). The data 
is saved in logarithmic size bins, each about 1.2 times 
larger than the previous one. The time is also measured 
logarithmically (next bin is 1.1 times larger than the pre- 
vious one). The extracted value for z in 3 dimensions is 
1.68 ± 0.07. The results for other dimensions are listed 
in Table | 
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H. Tables of Results 



measured exponents 


3d 


4d 


5d 


mean field 


1/u 


0.71 ±0.09 


1.12±0.11 


1.47 ±0.15 


2 


8 


0.015 ±0.015 


0.32 ±0.06 


1.03 ±0.10 


1 


(r + a08 - S)/au 


-2.90 ±0.16 


-3.20 ±0.24 


-2.95 ±0.13 


-3 


l/o 


4.2 ±0.3 


3.20 ±0.25 


2.35 ±0.25 


2 


t + a 08 


2.03 ±0.03 


2.07 ±0.03 


2.15 ±0.04 


9/4 


r 


1.60 ±0.06 


1.53 ±0.08 


1.48 ±0.10 


3/2 


d + 0/v 


3.07 ±0.30 


4.15 ±0.20 


5.1±0.4 


7 (at d c = 6) 


0/v 


0.025 ±0.020 


0.19 ±0.05 


0.37 ±0.08 


1 


avz 


0.57 ±0.03 


0.56 ±0.03 


0.545 ±0.025 


1/2 



TABLE I. Universal critical exponents. Values for the exponents extracted from scaling collapses in 3, 4, and 5 di- 
mensions. The mean field values are calculated analytically Jl2^|l3| ]. v is the correlation length exponent and is found from 
collapses of avalanche correlations, number of spanning avalanches, and moments of the avalanche size distribution data. The 
exponent 6 is a measure of the number of spanning avalanches and is obtained from collapses of that data, (r + a 08 — 3)/au 
is obtained from the second moments of the avalanche size distribution collapses. 1/a is associated with the cutoff in the 
power law distribution of avalanche sizes integrated over the field H , while r + a0& gives the slope of that distribution, r is 
obtained from the binned avalanche size distribution collapses. d + P/v is obtained from avalanche correlation collapses and 
0/v from magnetization discontinuity collapses, avz is the exponent combination for the time distribution of avalanche sizes 
and is extracted from that data. Error bars are based on variations in the results based on different approaches to the analysis: 
statistical fluctuations are typically smaller. 



calculated exponents 3d 4d 5d mean field 

a~08 0.43 ±0.07 0.54 ± 0.08 0.67±0.11 3/Z 

08 1.81 ±0.32 1.73 ±0.29 1.57 ±0.31 3/2 

0.035±0.028 0.169±0.048 0.252±0.060 1/2 

av 0.34 ±0.05 0.28 ± 0.04 0.29 ±0.04 1/4 

77 = 2 + (0-08)/v 0.73 ±0.28 0.25 ±0.38 0.06 ±0.51 



TABLE II. Values for exponents in 3, 4, and 5 dimensions that are not extracted directly from scaling collapses, but instead 
are derived from Table | and the exponent relations (see |l^|). The mean field values are obtained analytically |ll|[l3|. Both 
a0S and 08 could have larger systematic errors than the errors listed here. See the binned avalanche size distribution section 
for details. 



3d 4d 5d mean field 



Rc 


2.16 ±0.03 


4.10 ±0.02 


5.96 ±0.02 


0.79788456 


He 


1.435 ±0.004 


1.265 ±0.007 


1.175 ±0.004 





B 


0.39 ±0.08 


0.46 ±0.05 


0.23 ±0.08 






TABLE III. Non-universal scaling variables. Numerical values for the critical disorders and fields, and the rotation 
parameter B (equation in 3, 4, and 5 dimensions extracted from scaling collapses. The critical disorder is obtained from 
collapses of the spanning avalanches and the second moments of the avalanche size distribution. The critical field is obtained 
from the binned avalanche size distribution and the magnetization curves. H c is affected by finite sizes, and systematic errors 
could be larger than the ones listed here. The mean field values are calculated analytically [ p^|Jl^ ]. The rotation B is obtained 
from the dM/dH collapses. 
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V. COMPARISON WITH THE ANALYTICAL 
RESULTS 

Here we compare the simulation results with the renor- 
malization group analysis of the same system [ fb^lJTifl . Ac- 
cording to the renormalization group the upper critical 
dimension (UCD), at and above which the critical expo- 
nents are equal to the mean field values, is six. Close 
to the UCD, it is possible to do a 6 — e expansion, and 
obtain estimates for the critical exponents and the mag- 
netization scaling function, which can then be compared 
with our numerical results. 
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FIG. 14. Comparison between the critical exponents 
from the simulation and the e expansion Numerical val- 
ues (filled symbols) of the exponents r + g/38, t, 1/v, avz, 
and av (circles, diamond, triangles up, squares, and triangle 
left) in 2, 3, 4, and 5 dimensions. The empty symbols are 
values for these exponents in mean field (dimension 6). Ex- 
ponents in two dimensions are discussed elsewhere |p||l9|j2"o| . 
Note that the value of r in 2d conjectured value p5[~^ Wc 
have simulated sizes up to 30000 2 , 1000 3 , 80 4 , and 5(r7 where 
for 320 3 for example, more than 700 different random field 
configurations were measured. The long-dashed lines are the 
e expansions to first order for the exponents r + a0S, r, avz, 
and av. The short-dashed lines are Borel sums |34| for 1/v, 
as discussed in |p5[ . The lowest is the variable-pole Borel 
sum from LeGuillou et al. [Q, the middle uses the method 
of Vladimirov et al. to fifth order, and the upper uses the 
method of LeGuillou et al., but without the pole and with 
the correct fifth order term. The error bars denote system- 
atic errors in finding the exponents from extrapolation of the 
values obtained from collapses of curves at different disorders 
R. Statistical errors are smaller. 

Figure [l4| shows the numerical and analytical results 
for five of the critical exponents obtained in dimensions 
two to six (in six dimensions, the values are the mean 
field ones). The other exponents can be obtained from 
scaling relations jljj. The exponent values in figure [l4| 
are obtained by extrapolating the results of scaling col- 



lapses to either R R c or 1/L -> (see |j). The 
long-dashed lines are the e expansions to first order for 
t + a/38, t, avz, and av. The three short-dashed lines 
]i~3[ are Borel sums |Q for 1/v. Notice that the numer- 
ical values converge nicely to the mean field predictions, 
as the dimension approaches six, and that the agreement 
between the numerical values and the e expansion is quite 
impressive. 
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FIG. 15. Comparison between simulated dM/dH 
curves in 5 dimensions, and the dM/dH curve ob- 
tained from the e expansion. The thick dashed line shows 
the prediction of the e expansion to third order in e for the 
slope of the magnetization curve dM/dH in five dimensions. 
The theoretical curve is a parametric form j35| taken from the 
analysis of the ordinary, pure, thermal Ising model in three 
dimensions Il4] . The six simulation curves (thin lines) are for 
a system of 30 spins at disorders 7.0, 7.3, and 7.5 (R c = 5.96 
in 5 dimensions), and for a system of 50 5 spins at disorders 
6.3, 6.4, and 6.5. The latter curves are closer to the theoreti- 
cal dashed line). All the curves have been stretched/shrunk in 
the horizontal and vertical direction and shifted horizontally 
to lie on each other. 

The e expansion can be an even more powerful tool if 
it can predict the scaling functions. This has been done 
for the magnetization scaling function of the pure Ising 
model in 4 — e dimensions ]3^,|35| . Since the e expansion 
for our model is the same as the one for the equilibrium 



RFIM 1 13 1, and the latter has been mapped to all orders 
in e to the corresponding expansion of the regular Ising 
model in two lower dimensions fl3|,|37|,|38| , we can use the 
results obtained in 



35 1 . This is done in figure hf 



which shows the comparison between the dM/dH curves 
obtained in 5 dimensions and the predicted scaling func- 
tion for dM/dH, to third order in e, where e = 1 in 5 
dimensions (see [^5|). As we see, the agreement is very 
good in the scaling region (close to the peaks). 
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VI. SUMMARY 

We have used the zero temperature random field Ising 
model, with a Gaussian distribution of random fields, 
to model a random system that exhibits hysteresis. We 
found that the model has a transition in the shape of the 
hysteresis loop, and that the transition is critical. The 
tunable parameters are the amount of disorder R and 
the external magnetic field H. The transition is marked 
by the appearance of an infinite avalanche in the ther- 
modynamic system. Near the critical point, (Rc, He), 
the scaling region is quite large: the system can exhibit 
power law behavior for several decades, and still not be 
near the critical transition. This is important to keep in 
mind whenever experimental data are analyzed: decades 
of scaling need not imply self-organized criticality. 

We have extracted critical exponents for the magne- 
tization, the avalanche size distribution (integrated over 
the field and binned in the field), the moments of the 
avalanche size distribution, the avalanche correlation, the 
number of spanning avalanches, and the distribution of 
times for different avalanche sizes. These values are listed 
in Table ||, and were obtained as an average of the extrap- 
olation results (to R — ► R c or L — ► oo) from several mea- 
surements [l9[| . As shown earlier, the numerical results 
compare well with the e expansion [l^Jl^| . Comparisons 
to experimental Barkhauscn noise measurements p5| are 
very encouraging. 
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